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Abstract 



A two-replica graphical representation and associated cluster algorithm is 
described that is applicable to ferromagnetic Ising systems with arbitrary 
fields. Critical points are associated with the percolation threshold of the 
graphical representation. Results from numerical simulations of the Ising 
model in a staggered field are presented. For this case, the dynamic exponent 
for the algorithm is measured to be less than 0.5. 
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Monte Carlo simulations of equilibrium critical points have been revolutionized by cluster 
methods These methods effectively reduce critical slowing for a broad class of spin 

models. However, cluster methods developed thus far are efficient only for systems with a 
high degree of internal or translational symmetry. In this paper we show how to construct 
graphical representations and associated cluster algorithms appropriate for ferromagnetic 
Ising systems in the presence of fields. We demonstrate the method for the Ising system in 
a staggered field. 

The Swendsen-Wang (SW) cluster algorithm as applied to the Ising model updates 
both spin variables and bond variables. Starting from a spin configuration, satisfied bonds 
are occupied with probability p = 1 — e~^^ where (3 is the inverse temperature. A bond is 
satisfied if the spins at its two ends agree. Clusters of spins connected by occupied bonds 
are identified. Each cluster is randomly assigned a new spin value and each spin in the 
cluster takes that value. This constitutes one Monte Carlo step. The efficiency of the 
method is associated with the fact that the critical point of the Ising model coincides with 
the percolation point of the graphical model described by the bond variables Thus, at 
criticality, spin-clusters are coherently updated on all length-scales. 

A fundamental problem is encountered in applying the Swendsen-Wang method when 
fields are present. Clusters can be defined in the usual way but then the fields must accounted 
for in determining the probability of flipping clusters. Fields may be introduced directly via a 
Boltzmann factor or via ghost bonds. In either case, large clusters will tend to be 'frozen' 
in the sense that they almost always take one spin value. Because large clusters are frozen, 
equilibration occurs slowly on large length scales. In addition, the percolation transition of 
the graphical representation typically occurs in the disordered phase so that at criticality 
there are large clusters that are dynamically frozen. As a result, the qualitative gains 
expected from cluster methods (characterized by a small value of the dynamic exponent) 
are not realized. Nonetheless, limited quantitative success has been achieved for the random 
field Ising model by ad hoc methods that restrict the cluster size |Q. 

In this paper we present a graphical representation and associated cluster algorithm 
with the property that the percolation transition in the graphical representation coincides 
with the ordering transition in the spin-system. Furthermore, at criticality, the large-scale 
clusters are free to flip. We achieve this by using a replica method related to the replica 
Monte Carlo approaches that have been applied with some success to spin glasses p|-p!T|. 
We also note a related cluster method [p!^ -[I^, in which the system is "folded" on itself and 
pairs of sites are used to make clusters. 

The idea of our algorithm is as follows. Two independent replicas of the system in the 
same field are simulated simultaneously. Each site of the lattice is therefore associated with 

two spins, one from each replica, and can be in one of four spin states, (++), ( ), (H — ) and 

( — h). Sites where the replicas disagree, (H — ) and ( — h), are called active sites. Clusters of 
active sites are constructed and flipped. Allowed cluster flips interchange (H — ) with ( — h). 
Additional updating is applied independently to each replica to insure ergodicity. Since 
there is no net field on active clusters these flip freely. It turns out that percolation of active 
clusters signals the onset of long range order. 

The plan for the remainder of the paper is as follows. First we construct a joint distri- 
bution of the Edwards-Sokal [1^ type whose spin marginal is two independent Ising models. 
Next we indicate why percolation of the associated graphical representation coincides with 
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the onset of magnetic ordering in the spin-system. We then describe a cluster algorithm 
which simulates this joint distribution and present numerical results for the Ising model in 
a staggered field. Finally, we discuss generalizations of the method. 
The Hamiltonian for the Ising model is 

■^M = - ^'^i -^HiCJi (1) 

<i,j> i 

where the spin variables, cxj take the values ±1. The first summation is over the bonds 
of the lattice (or more generally, an arbitrary graph). The second summation is over the 
sites of the lattice and the fields are arbitrary. The Ising model on a square lattice in 
a staggered field (equivalently, the Ising antiferromagnet in a uniform field) is obtained by 
setting Hi = +H if i is in the even sublattice and Hi = —H if i is in the odd sublattice. 

We now define a joint distribution of two sets of Ising spin variables, {cTj} and {tj}, and 
a bond variable {rjij}. The bond variable is defined for each bond <i,j> and takes values 
and 1. We say that <i,j> is occupied if rjij = 1. The statistical weight X{a, r, rj) for the 
joint distribution is 

X{a,T,r]) = e"^^<-^>"""'"^"^+^^»'''^"'+"'^A(cr,r,r7)fip(r/). (2) 
B is the standard Bernoulli factor, 

BM=P^'KI-Pf'-^'^, (3) 

\ri\ = #{< i,j>\ r]ij = 1} is the number of occupied bonds and Nf, is the total number of 
bonds of the lattice. The A factor enforces the rule that only satisfied bonds are occupied: if 
for every bond <i,j> such that rjij = 1 the spins agree in both replicas (cr, = aj and Xj = Tj) 
then A{a,T,ri) = 1; otherwise A{a,T,ri) = 0. Equation is closely related to the "red- 
blue" graphical representation of the Ashkin- Teller model given in [l^ . It is straightforward 
to verify that integrating X{a, r, r/) over 7] yields the statistical weight for two independent 
Ising models in the same field, 

g-/3WH-/3HM ^ ^^^g^ V) (4) 

M 

if the identification is made that p = 1 — e~^^. 

Consider a two-replica spin system in which the a-replica has (-I-) boundary conditions 
and the r-replica has (— ) boundary conditions. The local order parameter is the difference 
between the magnetization of the two replicas, rrii = {< ai > — < Ti >)/2. Observe that 
magnetization in a single Ising model in a field is not generally the correct order parameter 
because the field induces local magnetization even in the disordered phase. By taking the 
difference between the magnetization of the two replicas with opposite boundary conditions, 
this contribution is canceled leaving only the spontaneous magnetization induced by the 
boundary conditions. 

Given a bond configuration rj we can ask for the conditional probabilities for the spins. 
Due to the A factor in the statistical weight, ai = +1 and Xj = — 1 if i is connected to the 
boundary by occupied bonds. On the other hand, due to the symmetry of the exponential 
factor in the statistical weight a site that is not connected to the boundary is equally likely 
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to be (H — ) or ( — h). Finally, (++) and ( ) spin states do not contribute to m^. Thus rrii 

is exactly equal to the probability that i is connected to the boundary and the onset of long 
range order coincides with percolation. For a more detailed argument, see |p!7 |. 



Our replica cluster algorithm simulates two independent Ising models, a and r, on the 
same lattice and in the same field. Sites i at which ctj 7^ Tj are called active sites. Bonds 
<i,j> connecting like spins in both replicas (cij = aj and r,, = tj) are called satisfied bonds. 

• Step 1: Satisfied bonds connecting active sites are occupied with probability p = 

1 - e-'f. 

• Step 2: Clusters of active sites connected by occupied bonds (including single active 
sites) are identified. The k^^ cluster is independently assigned a spin value, Sk = ±1 
with probability 1/2. If site i is in cluster k then the new spin values are ctj = Sk and 
Ti = —Sk- In this way all active sites are updated. 

• Step 3: Each replica is independently updated in a way that preserves detailed balance 
and insures ergodicity. This completes one Monte Carlo step. 

Without Step 3, the algorithm is not ergodic since the product (TjTj is locally conserved. 
Step 3 can be implemented in many ways. For example, each replica can be separately 
updated using the Metropolis algorithm. For the staggered field model in periodic boundary 
conditions we can effect further mixing by translating the r replica by a random amount 
relative to the a replica. If the translation is an odd vector, all r spins are fiipped. Since 
the Hamiltonian is invariant with respect to even translations and odd translations plus 
spin fiips it is clear that the translation part of the algorithm satisfies detailed balance. 
(Note that Metropolis sweeps are required even with random translations because the net 
staggered magnetization, s = [J2ieodd ~ J2i e even]{(^i + ^i) is otherwise a conserved quantity.) 
The simulations reported here implement Step 3 with both a Metropolis sweep and a random 
translation. 

The validity of the algorithm is proved by showing that it is ergodic and that the joint 
distribution, X{a,T,ri), defined in Eq. (||) is the stationary distribution of the algorithm. 



Ergodicity follows immediately from Step 3. To prove stationarity we observe, following |]T5 
that the Steps 1 and 2 of the algorithm correspond to conditional probabilities associated 
with X{a,T,T]). Step 1 is the conditional probability of a bond configuration given a spin 
configuration. Note that the bonds connecting inactive sites are not actually specified in 
the algorithm but since bonds are independently occupied this is of no consequence. Step 
2 is the conditional probability of the spin configuration on the active sites given a bond 
configuration, a set of active sites and the spin configuration on the inactive sites. 

We simulated the square lattice, staggered field Ising model in periodic boundary con- 
ditions using the two-replica cluster algorithm described above. Data was collected for 
f3H = 0, 2 and 4 and for size L in the range 16 to 256. Each L and (3H was simulated for 
50,000 Monte Carlo steps, dropping the first 5,000. Figure |l] shows the probability that there 
is a spanning cluster iS as a function of temperature T for several system sizes and (3H = 2. 
Spanning is defined as wrapping around the torus in either direction. The vertical line is the 



high precision value of Tc given in |T^. Figure ^ is a plot of S versus c{f3H)[T — Tc{f3H)]L 
of the data for all values of and L. Data collapse is achieved using Tc{(3H) from |1l8| . 
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c(0) = 1, c(2) = 2.64 and c(4) = 7.30. This figure illustrates that the model is in the Ising 
universality class independent of I3H. The fact that 5^1/2 for large systems and T < 
is due to periodic boundary conditions. For half of the Monte Carlo steps, replicas are 
magnetized in the same direction preventing active clusters from spanning. These results 
provide a clear numerical verification that percolation in the two-replica representation is 
coincident with the critical point. 

The Table shows the integrated autocorrelation time for the absolute value of the mag- 
netization of one replica, and the net staggered magnetization versus system size. The 
integrated autocorrelation time is 1/2 plus the sum of the normalized autocorrelation func- 
tion from time 1 through 200. The system size dependence of and Tm can be reasonably 
fit cither as AL^ or as A + B\og{L). For the whole range of L, logarithmic growth gives a 
better fit visually. Fitting a power law for system sizes greater than 40, we find dynamic 
exponents, Zm = 0.19 ± .09 and Zs = 0.33 ± .09 for 13H = 4 with nearly identical results 
for the other two field values. The quoted error is the statistical part and does not include 
systematic errors due to finite system size and finite cut-off in summing the autocorrelation 
function. A conservative conclusion is that O(log) < z < 0.5. It is clear that the algorithm 
achieves considerable acceleration over local dynamics, where z > 2. This is perhaps sur- 
prising in view of the fact that the algorithm uses local dynamic to break conservation of 
staggered magnetization. The autocorrelation times for the present algorithm and the ordi- 
nary Ising model {PH — 0) are roughly a factor of five larger than for the Swendsen-Wang 
algorithm however, further study is needed to determine whether the two algorithms share 
the same z. On the other hand, in the presence of a staggered field, autocorrelation times for 
both the Swendsen-Wang algorithm and the two-replica algorithm without translations are 
much larger than the values obtained here. Rough estimates of exponential autocorrelation 
times for the present algorithm show that they are about twice the corresponding integrated 
autocorrelation time. 

In the case of a staggered field, we have made use of the symmetries of the problem to 
incorporate as much mixing as possible into Step 3 of the algorithm. For systems such as 
the random field Ising model which do not enjoy translational symmetries, local dynamics 
are all that is available to equilibrate the average magnetization at each site. Thus, the 
two replica algorithm may not be a qualitative improvement over local dynamics alone. 
However, significant acceleration may be achieved by using many replicas. Suppose we have 
2K replicas, {a^^^\l = 1, . . . , K} all in the same field {Hi}. In each Monte Carlo step, the 
replicas are randomly paired and the two replica cluster procedure is applied to each pair. 
Each replica is also updated independently by some local ergodic algorithm. The replica 
summed magnetization at each site, Z^f^^ cf \ is conserved except by the local dynamics 
and so may relax slowly to equilibrium. This implies that the average magnetization at 
any site may reach equilibrium slowly resulting in a large exponential autocorrelation time. 
However, once the equilibrium values of the magnetization are reached, the fiuctuations of 
the replica summed magnetization are small for large K and thus couple weakly to the 
observables of individual replicas. This may yield rapid decorrelation and small values of 
integrated autocorrelation times. Note that in the two replica simulations, we observed Texp 
to be about twice Tint- 
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TABLES 



TABLE I. Integrated autocorrelation times for the absolute value of the magnetization of a 
single replica and the net staggered magnetization of both replicas Tg. For each entry, the one 
standard deviation error is 13%. 
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FIGURES 

FIG. 1. The spanning probability vs. T for various system sizes and (3H = 2. 

FIG. 2. The spanning probabiUty S vs. c{f3H)\T — Tc{(3H)]L for all system sizes and the three 
values of ^H. 
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